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Abstract: Image correlation spectroscopy (ICS) is known to be a useful 
tool for the evaluation of fiber width in the extracellular matrix. Here we 
evaluate a more general from of ICS fit parameters for fiber networks and 
arrive at a means of quantifying the fiber density, pore size and length 
which facilitates the characterization of the extracellular matrix. A 
simulation package was made to create images with different structural 
properties of fiber networks such as fiber orientation, width, fiber density 
and length. A pore finding algorithm was developed which determines the 
distribution of circular voids in the image. Collagen I hydrogels were 
prepared under different polymerization conditions for validation of our 
pore size algorithm with microscopy data. ICS parameters included 
amplitude, standard deviation and ellipticity and are shown to predict the 
structural properties of fiber networks in a quantitative manner. While the 
fiber width is related to the ICS sigma; the fiber density relates to the pore 
size distribution which correlates with the ICS amplitude in thresholded 
images. Fiber length is related to ICS ellipticity if the fibers have a 
preferred orientation. Findings from ICS and pore distribution algorithms 
were verified for both simulated and microscopy data. Based on these 
findings, we conclude that ICS can be used in the assessment of the 
extracellular matrix and the prediction of fiber orientation, width, density, 
length and matrix pore size. 

© 2012 Optical Society of America 

OCIS codes: (100.2960) Image analysis; (180.4315) Nonlinear microscopy. 
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1. Introduction 

Correlation spectroscopy uses autocorrelation to detect the likelihood of repetitions in a signal 
at a certain lag. In the case of two dimensional and spatially resolved data, this technique is 
called image correlation spectroscopy (ICS) [1,2] and has been used in numerous applications 
such as the estimation of fluorescent particle diffusion [3,4], degree of aggregation of plasma 
proteins [5] and the characterization of the extracellular matrix by analyzing the collagen fiber 
network. Previously Raub et al. [6,7] linked parameters from ICS analysis of multi-photon 
microscopy images with the bulk mechanical properties of collagen hydrogels. Here we 
evaluate how the properties of virtually computed fiber matrices affect parameters extracted 
from ICS and whether other parameters from the ICS analysis correlate with properties of 
fibrous networks. 

In ICS, once the autocorrelation function has been calculated, it is fitted to a suitable 
function to obtain a set of parameters. Here we discuss how the fit parameters of a more 
generalized Gaussian function correlate with variations in fiber density (fill factor), fiber 
length, fiber width and average orientation of the fibers. Our fit parameters include the 
amplitude, orientation, standard deviation (sigma), constant offset and ellipticity. With the 
help of simulated data, we establish a relationship between the fit parameters and the fiber 
properties. 

Raub et al. have shown how the collagen fiber width, length and pore size are related to 
polymerization conditions [7]. They found that ICS sigma correlates well with average fiber 
diameters while modest correlation was found in an attempt to estimate pore size based on 
ICS sigma [7]. The ICS's ability to create other parameters that are related to the voids within 
the network (pore size), fiber length and the average fiber orientation has not yet been shown. 

When employing fibrous networks in biologic experiments, fiber diameter and length are 
important parameters to characterize the network's mechanical properties. When evaluating 
the ability of cells or particles to migrate through the network, its pore size plays an important 
role [8,9]. In addition to the ICS technique, we have developed an algorithm through which 
we assess the pore distribution of the fibrous networks and we show that with appropriate 
image pre-processing steps the ICS fit parameters become correlated with the pore 
distribution. 

Our findings are further verified by applying the technique to second harmonic generation 
(SHG) images of collagen type I hydrogels. By using polymerization pH, we altered the 
structural properties of collagen gels [6,7,10] and applied our analysis technique to those 
images. 

2. Materials and methods 

2.1. ICS theory 

ICS entails the calculation of the autocorrelation function (ACF) of an image. The ACF is 
calculated more efficiently in the Fourier domain where the multiplication of the Fourier 
transform of the image (i) and its complex conjugate gives the power spectral density (PSD) 
and the inverse Fourier transform of the PSD the ACF as shown in Eq. (1): 

G(x, y) = T~ x \T{u{u,v)yT(i (11, v))] (1) 
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where u, v are pixel locations and x, y lags in pixels. As the ACF is concentrated about the 
center of the image, for structural analysis higher spatial lags can be removed by cropping a 
central portion of the ACF. Since the ACF at zero lag is related to noise, it is also removed in 
our analysis. In our case, we analyzed a central 12*12 portion of the ACF for a 256*256 pixel 
image. We determined the size of the 12*12 portion visually by comparing the ACF to a fitted 
function as illustrated in Fig. (1). The crop size should be minimized with the goal to cover 
large ACF values and representing the ACF drop-off with least amount of error. The crop size 
was kept constant for our study. 

2.2. Generalized Gaussian fit to describe the ACF 

Since the ACF is affected by the predominant orientation of the image's structural element, 
we modified the existing approach of parameter extraction [11] and included orientation as an 
additional analysis parameter. For illustration, Fig. (1) shows the dependence of the ACF on 
the oriented structural elements where those elements are fibers oriented primarily 
horizontally. In our modified approach, the ACF is fitted to a generalized 2-D Gaussian 
function shown in Eqs. (2). 

~ r \ a -( a ( x - x of+ 2b ( x - x o)(y-yo)+c(y-yo) 2 ) „ 

cos 2 9 sin 2 9 

a = — + — 

2<Tm 2a 2 m 

sin 20 sin 20 (2) 
b = — + — 

sin 2 9 cos 2 9 

c = — + — 

2^m 2a 2 




Fig. 1. Illustration of Image Correlation Spectroscopy and the dependence of the 
autocorrelation on orientation. The image with fibers oriented horizontally and its 
autocorrelation is shown. The Gaussian fit of the centrally cropped ACF is to the right with the 
wireframe representing the fit. 
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In Eqs. (2), x 0 , yo are the center of the ACF, a M , o m the standard deviation along the major and 
minor axis, theta the orientation of the axis with respect to the horizontal, A the peak 
amplitude and B the average offset. This allows us to summarize average structural properties 
of an image with five parameters. We term the ratio of the major axis standard deviation to 
the minor axis standard deviation as ellipticity which is another parameter in our ICS analysis. 

We modified an existing implementation of ICS [11] with the above function. The fitting 
is optimized in a least squares sense and the fit initialization conditions were adjusted to 
facilitate convergence. 

2.3. Simulation of images 

To verify the effectiveness of using ICS as a tool to provide a quantitative assessment of a 
fiber network, we simulated binary images with fibers of various properties. Lines represented 
fibers and were varied in length, thickness, orientation and density. Two types of random 
number generators were used, one creating values with a normal distribution (length and 
width) and the other with a uniform distribution (position, orientation). The starting point of 
each fiber was randomly selected over a range representing the desired image plus additional 
padding accommodating for average fiber length and variation. The orientation was uniformly 
selected in the specified range and an offset of 90 degrees was randomly added allowing the 
fibers to extend "forward" but also "backward" from the starting point. 

Table (1) describes the simulated parameter space. For each simulated condition we 
generated ten random images and results are reported as the average and standard deviation of 
those images. For example, we kept the mean fiber length, fiber density constant and varied 
the mean fiber width with the fibers following a normal distribution. This gave us images to 
evaluate the influence of fiber width on the fit parameters. This approach was repeated for 
varying mean fiber content and length. To evaluate fiber length, fibers oriented with different 
uniform distributions ranging from a dominant to completely random orientation were 
simulated for varying fiber lengths. This enabled us to show the dependence of fiber 
orientation and length on ellipticity. 

Table 1. Simulation variables 
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2.4. Method for estimating the pore distribution 



Since our fiber network simulation program produced a random image without knowing the 
void size and distribution we developed a technique to estimate the average void size. The 
void size is approximated by finding the largest circular object fitting into a void (pore). 

The flowchart to calculate the pore distribution is shown in the Fig. (2a) and can be 
illustrated by the following example: The input image is a binary image which has been 
simulated with a fill factor of 0.5, mean fiber length of 50 pixels, mean fiber width of 2 pixels 
and the orientation following a uniform distribution between [-90, 90] as shown in Fig. (2b). 
A circular mask with a diameter of 10 pixels is convolved with the image (Fig. 2b). Areas 
with zero intensity are identified in the convolved image, which corresponds to voids where 
the mask is capable of fitting. Figure (2c) illustrates all regions with intensity of zero (black 
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area) however, independent regions with an area of one pixel are the locations where the 
circular mask fits exactly (marked with crosses). When increasing the mask size, the 
independent regions shrink. They often do not shrink to exactly one pixel and the smallest 
area is several pixels in size. Therefore we choose an arbitrary threshold of 10 pixels and 
recorded the centers of all independent regions with a size smaller than this threshold. 
Independent regions were identified with connected component labeling (Image processing 
toolbox, Matlab, MathWorks, Natick, MA) [12] which provides the centroid and area for each 
region. We carried out the entire process for masks representing all diameters of interest and 
maintained a list of possible pore centers. Centroid coordinates were rounded to the nearest 
integer pixel and if more than one pore size was found at the same location, the one associated 
with the largest diameter was kept (Fig. 2d). The list of centroids and diameters enables us to 
create the pore distribution as shown in Fig. 2(e). Our experimental results indicated that a 
lognormal function fits the pore distribution best. 




Fig. 2. Flow chart and illustration of pore finding algorithm. (A) Steps entailed in finding the 
pore distribution of an image. The algorithm produces a list of pore centers and the associated 
diameters. (B) The image is convolved with a circular mask. To find the pore distribution, 
circular masks with different diameters are generated. Zeros of the convolved image are 
identified which corresponds to the points where the circular mask fits into the pore. (C) 
Connected component labeling is used to identify independent regions that have an area 
smaller than a set threshold. For illustrative purpose, centers of those regions are marked with 
crosses (C left) and the corresponding circular pores are drawn in the subfigure to the right. (D) 
Simulated image with all pores identified. (E) Pore diameters [pixel] follow a lognormal 
distribution. 

2.5. Microscopy of collagen type Ihydrogels 

Imaging was performed with a single beam intravital microscope (TrimScopel, 
LaVisionBioTec, Bielefeld, Germany). The Titanium Sapphire laser was tuned to 780nm and 
the field of view was about 350 microns. SHG images of collagen I hydrogels were recorded 
similar to previous publication of our group [13]. An image stack was formed by varying 
focus in the sample and five consecutive images were analyzed for each sample. The average 
and standard deviation of those analyses is reported for each analyzed hydrogel. 5ml collagen 
mixture was prepared and with 35ul lMNaOH added vortexed to mix. 250ul mixture was 
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removed and added to the first well of a 48 well plate sitting on ice. 5ul NaOH was added to 
the remaining mixture and vortexed, while another 250ul was removed to the second well. 
This was repeated for 13 wells and the volume of NaOH was reduced to 2.5ul for wells 14 
through 15 resulting in collagen gel pH varying from 6 to 9. For this study repeatable 
polymerization outcome under different pH conditions was not necessary as we primarily 
wanted to study hydrogels with varying structural properties. 

2.6. Preprocessing of microscope images for ICS and pore finding algorithm 

Our microscopy grayscale data was subjected to Wiener filtering for noise removal (3x3 
kernel, Image processing toolbox [14], ) which improves the output of the pore distribution 
algorithm for our imaging data. The raw grayscale images are not suitable for thresholding 
directly as the intensity on the optical center is brighter than at its edges. Adaptive histogram 
equalization (Image processing toolbox) [15] was applied to create uniform grayscale 
distribution in 8x8 sub areas of the image (Fig. 3). Alternatively the images could be corrected 
for nonuniform illumination by measuring the spatial system response. In order to create a 
binary image, the images were then thresholded using Otsu's method (Image processing 
toolbox) [16]. As shown in the flowchart in Fig. (3), once the binary image is obtained we 
conduct ICS and Pore distribution analysis. 

It is possible to conduct ICS directly on a gray scale image. We have tested ICS on raw 
data, adaptive histogram corrected gray scale images and on threshold images and found that 
for comparing ICS with our pore routine, analysis on thresholded and histogram adjusted 
images created most reproducible results (data not shown). 



Image after adjusting the contrast Adaptive histogram equalization 




Fig. 3. Steps involved in ICS and pore distribution analysis of measured data. The pore finding 
routine and ICS is applied on images after contrast adjustment, adaptive histogram equalization 
and thresholding. Image (1) illustrates the raw data, (2) the histogram adjusted data, (3) the 
thresholded image, and (4) illustrates the pore centers and diameters. 



3. Results 

Representative images simulated for different widths, fiber densities and lengths are shown in 
the Fig. (4). ICS sensitivity analysis was conducted on those images as well as the pore 
finding algorithm was applied to identify if any of the ICS parameters correlated with the pore 
diameter. 
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Width = 1 Width = 5 +/- 1 Width = 9 +/- 1 




Fill factor = 1 0% Fill factor = 50% Fill factor = 80% 




Length = 5 +/- 1 Length = 50 +/- 5 Length = 200 +/- 20 




Fig. 4. Images of virtually computed fiber matrices under different simulation conditions. The 
top row contains images of different widths with a constant mean fiber density and length. The 
middle row has images of different fiber densities with constant mean width and length. The 
bottom row illustrates images of different lengths with constant mean fiber density and width. 



3.1. ICS parameter sensitivity 

When plotting ICS fit parameters versus fiber properties, it can be observed in Fig. (5) that 
with increase in fill factor; the amplitude, sigma and offset decreases. However the change in 
amplitude was more strongly affected by fill factor then the other parameters. Similarly an 
increase in fiber width increases amplitude, sigma and offset; however the increase in sigma is 
more pronounced than amplitude and offset. It can also be observed that while the fiber length 
correlates weakly with the amplitude and offset, it is a strong function of ellipticity when 
fibers show a preferred orientation. We can conclude that the amplitude is influenced the most 
by the fill factor, sigma by the width and ellipticity by the length of the fibers subject to a 
dominant orientation. When comparing amplitude, sigma, offset and ellipticity versus 
preferred orientation, we found that amplitude is slightly increased with preferred orientation 
(data not shown), but the change is more than an order of magnitude smaller than the change 
induced by the fill factor. Offset was very little affected by preferred orientation as changes 
were within the error of our simulations. Ellipticity was strongly affected by preferred 
orientation because sigma along the major axis was strongly affected by preferred orientation. 

3.2. Pore size distribution 

The simplest approach to simulate a fiber network with smaller pore size is to increase the fill 
factor in our simulations. We calculated the pore size distribution for different fill factors 
ranging from 10% to 80%. We found that the fill factor correlates well with average pore size 
which was obtained from a log normal fit as shown in Fig. (2e). We also observed that the 
pore size variation within an image did not show a trend with the fill factor (data not shown). 
Then ICS analysis was conducted on the same images to confirm whether ICS parameters 
correlate with the pore size. We found that among all the ICS parameters, ICS amplitude 
varies the most with changing mean pore diameter (Fig. 6, left). That sensitivity was stronger 
than the influence of offset and any of the other ICS parameters (data not shown). 
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Fig. 5. ICS parameter sensitivity analysis. The top left figure illustrates sensitivity of fit 
amplitude with fill factor, width and length. Symbols represent and average of 10 simulations. 
Standard deviation is also plotted but often smaller than the symbol labeling the plot. The 
amplitude is highly correlated with fill factor. The top right figure illustrates the variation of 
average sigma, with sigma shown to be highly correlated with fiber width. The average sigma 
is the mean of sigma on the minor and major axis. The bottom left figure illustrates the offset 
with offset shown to be highly correlated with fill factor. The bottom right figure illustrates 
ellipticity to be highly correlated with fiber length if the fibers have a preferred orientation. 
Orientations follow uniform distributions from [-5°, 5°] (Orient = 5) to [-80°, 80°] (Orient = 
80). Uniform distribution of [-80°, 80°] indicates almost randomly oriented fibers and the 
ellipticity (o M /o" m ) is almost flat around 1 . 




Pore Size [pixel] Pore Size (e^ f ) [pixel] 

Fig. 6. Correlation of ICS amplitude with average pore size of simulated (left) and microscopy 
data (right). The microscopy data (n = 15) spans a pore size range that is much smaller than the 
simulated data. 
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3.3. ICS parameters from microscopy data 

In order to verify our findings on microscopy data, we prepared collagen gels under different 
polymerization conditions. Our raw grayscale images showed that for low and high pH 
values, the SHG signal level is low and higher signal levels were obtained when 
polymerization was carried out at room temperature as compared to 37 degrees C (data not 
shown). The raw grayscale images were preprocessed as described above and subjected to 
ICS as well as pore distribution analysis. After preprocessing the fill factor varied between 0.5 
and 0.8. For that limited range, the ICS amplitude correlated as expected with fill (R 2 = 0.92) 
and ICS offset had little correlation with fill factor or (R 2 = 0.23) and sigma (R 2 = 0.36) 
although trends were consistent with the simulated data (data not shown). ICS amplitude was 
plotted versus average pore size (Fig. 6, right) and it can be observed that the mean pore 
diameter follows a linear trend with the ICS amplitude. Also, our experimental results 
confirmed that the pore size had the best correlation with ICS amplitude as previously found 
in the simulated data. Using a linear model, the goodness of fit (R 2 = 0.901) indicated that the 
pore size correlates reasonably well with the ICS amplitude as opposed to for example the 
ICS offset (R 2 = 0.01). 

4. Discussion and conclusion 

With the results on simulated images we have demonstrated how our ICS parameters correlate 
with the structural properties of virtual fiber networks. While the Gaussian fit was dependent 
on the orientation of the structural elements in the image, the width of the fibers was found to 
strongly correlate with the sigma. The fiber length is a major structural parameter and we 
found a strong dependence of fiber length on ellipticity when fibers tend to align in a 
particular direction. In simulated data amplitude is decreasing the most at small fill factors 
while ellipticity changes the most for shorter fibers. 

When comparing the ICS fit parameters from microscopy data we found the same trend as 
observed in our simulated data. There was a strong correlation of ICS amplitude with fill 
factor. We did not find correlation between ICS offset and any of our computed parameters. 
Our microscopy data covered a fill factor of 0.5-0.8 where ICS amplitude change is smaller, 
never the less we were able to obtain a high linear correlation coefficient. Ellipticity in our 
microscopy data varied between 1 and 1.6 indicating we had a small preferred orientation 
with relatively short fibers. In that range ellipticity changes the most which might explain why 
we could observe an increase in ellipticity with increased pH above pH 7 (data not shown). 
Others have reported that with increase in the basic environment, fibers become longer [7] 
which can lead to an increased ellipticity if there is a preferred orientation in the field of 
view.. This is interesting as a preferred orientation is not obvious in our images. This suggests 
that only a small deviation from uniform orientation is necessary to observe this effect. 
During cell migration in collagen hydrogels, it is likely that cells will create preferentially 
oriented fibers along their migration track and ICS ellipticity might be useful in analyzing 
fibers along those tracks [13]. Additionally, collagen fibers in tissue are generally not 
randomly oriented and ICS ellipticity might be a useful quantitative parameter for tissue 
imaging. 

Our pore finding algorithm is an effective method to describe pore size distribution. 
Because ICS amplitude decreases with an increase in fiber density we were able to compare 
the pore size results with ICS. We found a correlation between the mean pore size and the ICS 
amplitude. This has been validated with our simulations as well as our SHG image data. Pore 
size of our microscopy data was at the lower range of the simulated data where it is more 
difficult to observe a change, however we still found a good correlation coefficient. In order 
for the ICS amplitude to correlate with pore size it was important to correct the images for 
nonuniform system response and to threshold our data. 

All our analysis was conducted on two dimensional images. The ACF and the Gaussian fit 
can be expanded to three dimensions. Similarly the pore finding algorithm can be expanded 
with a third dimension. While analysis with a third dimension will not change the ability to 
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assess fiber diameter, it would affect predicted fiber length. In confocal microscopy the 
resolution along the optical axis is about three times worse than the planar resolution which 
allows estimating fibers' lengths even when they are at oblique angles in respect to the 
measurement plane, however their length is underestimated. When comparing microscopy 
images in Fig. 3 with the simulated images in Fig. 4 one can observe that fibers can exhibit 
tapering towards the end as well as splaying. Our line algorithm does not simulate these 
characteristics but it is likely that they lead to an underestimation of fiber length. 

We can conclude that the structural properties of collagen fibers have been evaluated with 
Image Correlation Spectroscopy and we have shown that ICS parameters correlate with fiber 
width, length, average orientation and pore size of simulated fiber networks. Pore size 
calculations were verified on collagen I hydrogels and showed good correlation. 
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